function sa2c(N)

  

  
  time_steps = [0:N-1]./(N-1)*0.999;

  function resSa2c = Sa2c(z)
    if (z < 0.001)
      resSa = 1+(-9/8).*z+(5/4).*z.^3+(-9/7).*z.^4+(4/3).*z.^6+(-27/20).*z.^7+(11/8).*z.^9+(-18/13).*z.^10+(7/5).*z.^12+(-45/32).*z.^13+(17/12).*z.^15+(-27/19).*z.^16+(10/7).*z.^18+(-63/44).*z.^19;
    else
      resSa = z.^(-2).*((-1/2)+(1/4).*(2.*((-2)+(-2).*z+z.^2).*(1+z+z.^2).^(-1)+z.^(-1).*(3.^(1/2).*pi+(-2).*3.^(1/2).*atan(3.^(-1/2).*z.^(-1).*(2+z))+3.*log(1+z+z.^2))));
    end


  end
  function resM2c = M2c(z)
  % pulled out z^3
    resM2c = zeros(2,2);
    resM2c(1,1) = 1;
    %resM(2,1) = 2 + z.^3;
    %resM(2,2) = z.*(z.^3 - 1);
    resM2c(2,1) = (-4+7.0.*z^3);
    resM2c(2,2) = z.*(z.^3-1);
  end
  function res2c = deq2c(z, u)
    %disp(u)
    a   = u(1);
    az  = u(2);
    res2c = zeros(2,1);
    res2c(1) = az;
    res2c(2) = 2*Sa(z) - 9.0.*z.^2.*a - 2.0.*z;
  end
  %uin = [(2+pi./3./sqrt(3)-log(3))./3; -4./6];
  uin = [(-1/3)+(1/36).*(18+3.^(1/2).*pi+(-9).*log(3));1./4];
  uin
  opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-10, 'MStateDependence', 'none', 'Mass', @M2c);
  [tout, uout] = ode15s(@deq, time_steps, uin, opts);
  %plot(tout, uout(:,1)./log(1-tout))
  plot(uout)

  %uout(:,1)./log(1-tout)


end
